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Abstract - The computational fluid dynamics code FLUENT has been 
adopted to simulate the entire rectangular-channel-like (3-D) geometry 
of an experimental CVD reactor designed for Si deposition. The code 
incorporated the effects of both homogeneous (gas phase) and 
heterogeneous (surface) chemistry with finite reaction rates of 
important species existing in silane dissociation. The experiments were 
designed to elucidate the effects of gravitationally-induced buoyancy- 
driven convection flows on the quality of the grown Si films. This goal 
is accomplished by contrasting the results obtained from a carrier gas 
mixture ..of .I^/Ar with the ones.: obtained from* .the. same molar mixture 
ratio of ^/He, without 'any accompanying change in“ the chemistry. 
Computationally , these ’cases are simulated in 'the terrestrial, 
gravitational field and in the. absence of gravity. The numerical 
results compare -favorably with experiments'. ‘Powerful computational 
tools provide invaluable insights ‘into the complex physicochemical 
phenomena taking place in CVD reactors. Such information is essential 
for the improved design and optimization of future CVD reactors. 

1 - INTRODUCTION 

Despite its routine application for the production of many electronic, optical 
and structural materials, the level of understanding in CVD processes has been 
rather limited. This is evident from the conventional approaches, taken by 
CVD reactor designers, which are mostly based on intuition and experience. 

The sometimes unjustified assumptions and simplifications employed by many 
modellers are also consequences of insufficient fundamental knowledge about 
the. complex physicochemical interactions involved. Extensive discussions of 
available CVD models are given in the recent reviews of Hess, Jensen and 
Anderson/1/ and Jensen/2/ . Indeed, a comprehensive analysis should 
simultaneously include, among other things, homogeneous ( gas phase) and 


*NASA Senior Resident Research Associate. 


i 



heterogeneous (surface) chemical reaction kinetics, heat and multicomponent 
mass transport, fluid physics and thermodynamics. In many cases, these 
phenomena should be treated in multidimensions, and, in some, even the 2-D 
approaches fail to provide the necessary insight, as is demonstrated by 
Rosenberger/3/ . 

There are, however, some reactor/flow/susceptor configurations and thermal 
conditions where certain approximations can be made. For such cases, an extra 
level of sophistication is not warranted. Stagnation flow reactors studied by 
Rebenne and Pollard/4/, Houfman, Graves and Jensen/5/ and Wahl,et al./6/, or 
rotating disk reactors studied by Evans and Grief /7/ and Coltrin et al./8/, or 
horizontal reactors studied by Coltrin, Kee and Miller/9,10/ and Ristorcelli 
and Mahajan/11/ are examples where the analyses are justifiably reduced to 2-D 
with negligible sacrifice in information. 

Recently, increasing interest in the "CAD-CAM" approach to reactor design 
methodologies and improving availability and quality of homogeneous chemical 
kinetic data have to motivated full 3-D transport analyses. The application 
of the computational package PHOENICS by Rhee, Szekely and Ilegbusi/12/ to the 
analysis of a rectangular reactor, neglecting both homogeneous and 
heterogeneous finite rate chemistry, and the systematic study of transport 
phenomena in horizontal and tilted CVD reactors by Moffat and Jensen/ 13/ 
demonstrate such trends. 

Presently, the computational fluid dynamics code FLUENT is being enhanced and 
extended under a NASA program to incorporate many of the phenomena encountered 
in CVD. The intermediate results of this development effort are reported in 
references 14 and 15. We have adopted this intermediate version of FLUENT-CVD 
to study the effects of both buoyancy and chemistry in the Si deposition 
reactor designed by Tsui/16/ (see also Tsui, P. and Spear, K.E., to be 
published) . 

The nearly two-decade old Si deposition rate results of Eversteyn et al./17/ 
have so far been the only experimental data for modellers to test their codes 
against. Relevant and reliable CVD rate data is hardly available in the 
literature for the validation of sophisticated models under stringent 
conditions. Tsui* s/16/ Si deposition experiments were, however, designed to 
address a modeller's needs. Therefore, they provide excellent benchmark cases 
for code verification. 

Tsui's experiments were designed to elucidate the the role of 
gravitationally-induced buoyancy-driven convection flows on the CVD rates and 
the quality of the Si films grown. Thus, rates obtained from different 
carrier gas mixtures of H 2 /He are compared with the ones obtained from 
corresponding molar mixtures of H 2 /Ar. This assures that no accompanying 
change in the chemistry is introduced. The experiments were performed at 
various substrate temperatures using 0.1 mole percent silane as the Si-source 
gas. Previous modelling efforts of Stinespring and Annen/18/ for these 
experiments predicted qualitative and quantitative trends which did not agree 
with the experimental observations. This can be mainly attributed to a) their 
model assumes fully developed velocity and temperature fields at the hot 
substrate leading edge with no entrance length for their development, b) their 
species equations are only 2-D in a 3-D flow field, c) their computational 
volume takes the side edges of the susceptor as the location of reactor side 
wall which oversimplifies the real situation (indeed, this assumption is 
routinely employed by even the 3-D modellers in the literature with sometimes 
serious consequences, i.e. different temperature fields, aspect ratios, 
buoyancy effects, etc., and is not a realistic description of typical CVD 
reactors), and d) they assume thermochemical equilibrium (infinitely fast 
reaction rates) in the gas phase. 

The main objective of this study is to address the question of how significant 
the effect of buoyancy is on the quality of the crystals grown by CVD. Our 
preliminary comparison with experiments indicate that FLUENT-CVD can correctly 
describe many of the 3-D physicochemical phenomena occur ing in a typical CVD 
reactor. However, more rigorous tests are needed to confidently validate the 
code. The results of such a thorough comparison will be reported elsewhere 
(Gokoglu et al. , to be published). 
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2 - PHYSICAL MODELS AND MATHEMATICAL FORMULATION 


The details of the transport equations used by FLUENT- CVD describing the 
conservation of mass, momentum ( Navier-Stokes ) , energy and species mass are 
given in references 19 and 20. In Cartesian tensor notation they can be 


written as : 

Total Mass : 

a/axjjpu-L) =o (1) 

Momentum : 

3/3x t (pu i u a ) = -3p/3x ± + S/SxJjiOuj/SXfc+SUit/SXi) ] + pg ± (2) 

Energy : 

a/axjjpu-jh) = 3/3xj_ (A/Cp’3h/ 3x^ ) (3) 

Species Mass : 

S/SXjJpUjYj) = 3/3Xj_(p/Sc j *3Yj/3x^) + 1^=1 Rj^ (4) 

The temperature, T, is obtained from the enthalpy, h, via the relationship 
dh = CpdT. The mixture density, p, is calculated from the ideal gas law using 
the equation : 

P = p/fRTZjLi Yj/Wj) (5) 


where J is the total number of gaseous species. The dynamic viscosity, p, and 
the thermal conductivity. A, of the mixture are obtained from the 
semiempirical formula of Wilke/21/. The mixture specific heat, c p , is 
calculated from the molar weighted sum of the individual species specific 
heats. Polynomial curve fits in T are subsequently calculated to describe the 
temperature dependences for each of p, A and Cp to be eventually used by 
FLUENT-CVD. Natural convection is included via the gravitational body force 
in Eq . ( 2 ) and the temperature dependence of p given by Eq.(5). The Schmidt 
number, p/(pDj), is assumed to be constant for each species. Therefore, the 
diffusion coefficient, D j , will have the same temperature dependence as p/p. 

It should be noted that currently Eq.(4) does not include thermal (Soret) 
diffusion of species. This is certainly not a justifiable assumption for most 
CVD applications where steep temperature gradients exist/22,23/. Further 
development work of FLUENT-CVD includes plans to incorporate thermal diffusion 
in the near future. 

It should also be noted that because the Si-carrier species are dilute, 

( Yj<< 1 ),in the host gas mixture, a number of simplifications can be made, 
sucn as 1) neglecting the effect of Si-carrier species compositions on fluid 
properties, p, A and c p , and the species diffusion coefficient, Dj , 2) using 
the binary, rather than the multicomponent, diffusion coefficients for each 
species in the mixture , 3 ) neglecting the rate of creation of energy by 
homogeneous chemical reactions in Eq. (3) and by surface reactions, and 4) 
neglecting the energy transport associated with species mass fluxes (Dufour 
effect ) in Eq. ( 3 ) . 

Homogeneous Reactions : 

Rather than employing the complex system of reactions for Si deposition from 
silane dissociation involving more than 20 different molecules, as reported by 
Coltrin, Kee and Miller/9,10/, we have limited the significant gas phase 
molecules to silane, silylene and disilane, based on the discussion given by 
Moffat and Jensen/13/. Indeed, since there is no experimental evidence to 
suggest otherwise, this approach is justified in hydrogen rich environments, 
as is the case studied here. Therefore, the reduced reactions are given by : 
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SiH 4 


k l,f 

k l,r 


SiH 2 + H 2 


(Rl) 


k 2 f 

Si 2 H 6 =' SiH 4 + SiH 2 
k 2,r 

An Arrhenius type rate expression, k=ATPexp( -E/RT) , 
and reverse reaction rate constants, given by/13/ 


(R2 ) 


is used for the forward 


Reaction 

Molar 

Mixture 

A(s 1 ) 

A(m 3 /kg-s) 

3 

E( J/kgmol) 

i,f 

all 

6.10-10 28 



-5.00 

2 . 4614 • 10 8 

2 , f 

all 

2.12-10 35 

-- 

-6.47 

2.3597 • 10 8 

l,r 

25%He/75%H 2 

— 

2.087-10 21 

-4.44 

1.4267-10 7 

l,r 

25%Ar/75%H 2 

— 

4.583-10 20 

-4.44 

1.4267-10 7 

2,r 

25%He/75%H 2 


7.075- 10 23 

-4.50 

1.2845-10 7 

2 , r 

25%Ar/75%H 2 

— 

1.554-10 23 

-4.50 

1 . 2845 • 10 7 


Modelled Reactor and Boundary Conditions : 

The modelled portion of the schematic overview of Tsui 1 s experimental reactor 
is given in Figure 1. Full details of the actual reactor are given in 
reference 16, The susceptor is made of silicon carbide coated graphite and is 
inductively heated. The substrate is a 5.08cm long by 2.54cm wide silicon 
wafer centered on the susceptor. The solid block on which the susceptor is 
placed and the block upstream of the susceptor are both made of quartz. These 
solid bodies are naturally heated by heat conduction from the hot susceptor. 

A surface temperature of 1373K is used for the heated susceptor top, side and 
back surfaces. The reactor unheated (cold) surfaces are assumed to be at 300K 
everywhere except at the bottom quartz surface of the rectangular entrance 
duct upstream of the heated susceptor. The conjugate heat transfer capabilty 
of FLUENT-CVD is utilized in that region to actually calculate the steady 
state equilibrium surface temperatures resulting from heat conduction along 
the bottom surface of the entrance region prior to the deposition surface. 

In order to simulate the experimental conditions, we have used the same 4 slm 
total carrier gas (75% H 2 with 25% He or 75% H 2 with 25% Ar by mole) flowrate 
with 0.1 mole percent silane at 300K for the reactor inlet conditions. At the 
reactor exit, the program considers global mass conservation and imposes a 
zero normal gradient condition for the solution variables. On the solid 
boundaries, the non-slip condition is used for velocity components parallel to 
the walls. Because the Si-carrier reactive species are dilute in the mixture, 
a zero convective flux condition is used normal to the walls. 

As the boundary condition for species on heated surfaces we use the balance 
between the diffusional mass flux of each species to/from the surface and the 
local rate of creation/destruction of that species due to chemical reactions 
at the surface. Thus, for each chemical species, j, 

K s 

“pDj n*V Yj - Rjk (6) 

at the surface. Note again the absence of the thermal ( Soret ) diffusion term 
in the species diffusional mass flux expression. Although Eq.(6) applies on 
all surfaces, for practical computational efficiency purposes a zero species 
mass flux condition is used on the unheated walls with negligable loss in the 
accuracy of results. The growth rates of Si on heated surfaces can then be 
obtained from the deposition rates by 

-j K s 

Growth rate(length/time) = [p S j_( s ) Rsjjsjk (7) 
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The jth surface reaction rate can be expressed as the product of the molecular 
flux across the Knudsen sublayer (estimated from the frequency of molecular 
surface collisions using the kinetic theory) and the reaction probability, ^ , 
(sometimes referred to as the sticking coefficient) of the jth species J 

Rj = 'rj[RT/(2 1 rW j )] 1 /2 p Y j (8) 

Because FLUENT-CVD requires the surface reactions expressed in terms of 
Arrhenius type [ AT^exp ( -E/RT) ] rate constants (unless they are provided by 
user supplied subroutines), Eq.(8) is rewritten in the form of a product of a 
pseudo-first order reaction rate constant and concentration (pYj) for each 
species. We have also used the surface reaction probability expression of 
Coltrin, Kee and Miller/10/ for silane given by 

^SiH4 =5 * 37 * 10“ 2 exp( -9400/T) (9) 

which can be combined with the Arrhenius rate constant expression. The 
surface reaction probabilities of the other reactive Si-carrier species, SiH 2 
and S^Hg , are taken to be unity whereas those of H? and the inert gases to be 
zero. The resultant A, £ and E values are tabulated below. 


Species 

Atm/s-K 1 /^ ) 

P(K 1/2 ) 

E ( J/kgmol ) 

SiH 4 

0.3453 

0.50 

7.8156* 10 7 

SiH 2 

6.6412 

0.50 

0.0 

Si 2 h 6 

4.6197 

0.50 

0.0 

*2 

0.0 

- 

- 

He 

0.0 

- 

- 

Ar 

0.0 

- 

- 


3 - RESULTS AND DISCUSSION 

The deposition of silicon takes place on any surface where the temperature is 
high enough for surface reactions to occur. Therefore, a fair comparison of 
Si deposition rates to experiments requires that the model keeps a correct 
account of the Si-carrier species inventory throughout the reactor, contrary 
to what most models in the literature do, i.e. to consider the susceptor or 
the substrate region as the only deposition surface. Thus, when calculating 
the Si deposition rates on the susceptor top surface, the code considers the 
deposition on all relevant surfaces, including the quartz surface upstream of 
the susceptor and the vertical side and back surfaces of the susceptor. 
However, the top substrate surface is the focus of this study for comparison 
with the experiments. Therefore, the presentation of the results will mostly 
concentrate on the susceptor area with the three X-slices (Y-Z plane) shown in 
Figure 2 at the leading edge, near the center and at the trailing edge of the 
susceptor. 

Because of the Z-symmetry of the reactor assembly, the computational domain 
uses only a half of the reactor and exploits reflection boundary conditions on 
the central Y-Z plane. For the finite-difference formulation of FLUENT-CVD, 
we used a computational mesh for the half-reactor which had a total of 31 
cells in X, 16 cells in Y and 10 cells in Z directions. The grid is 
nonuniformly distributed to accomodate regions of large gradients in solution 
variables, especially immediately above the susceptor. Even with this 
computational scheme, the solutions presented in this paper are still grid- 
sensitive. However, it was ascertained that the correct physicochemical and 
transport phenomena are resolved. A finer mesh is not warranted with the 
current limitations of FLUENT-CVD and final fully converged results will be 
published at a later time. 

This study uses the silane dissociation chemistry for silicon deposition. A 
25 mole percent of either Ar (referred to as the Ar case hereafter) or He 
(referred to as the He case hereafter) in 75 mole percent H 2 is used as the 
carrier gas mixture at a constant volumetric flowrate of 4 slm. These two 
distinct cases were chosen by Tsui because of the presumed strong buoyancy- 
driven convection for Ar and the much reduced level of such convection for He. 
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The added silane (0.1 mole percent) replaces hydrogen. Note that because the 
mole fractions are fixed and inert gases replace hydrogen, no accompanying 
changes in the chemical reaction mechanisms are expected. 

Considering the couplings among the forced-convective, free-convective and 
diffusive transport processes in chemically reacting flow systems, choosing a 
criterion by which unambigious information can be extracted for the effect of 
only free-convection from terrestrial experiments is not trivial. Depending 
on the focus of interest, different experimental design conditions ( such as 
iso-f orced-convection , i.e. constant Reynolds number, iso-forced-convection- 
diffusion,- i.e.‘ constant species surface flux or deposition rate when gravity 
is absent, iso-mass flowrate, iso-volumetric flowrate, etc.) can be employed 
as the basis for comparison. It should be pointed out that Tsui's experiments 
were run at constant total volumetric flowrates. This means that the forced- 
convection, conditions (Reynolds numbers) and species dif fusivities (Schmidt 
numbers) were .different for the two Ar and He cases. Indeed, under the 
conditions of Tsui's experiments, neglecting the surface kinetic barriers, the 
He case is estimated to give about 13 percent higher deposition rates than the 
Ar case even' when the buoyancy effects are absent. Therefore, the differences 
between the terrestrial experimental deposition rates for the Ar and He cases 
should be studied accordingly, i.e. if the Si deposition rates (and their 
distribution over the substrate) are chosen to be the experimentally 
observable-parameter to base the conclusions on, the interpretations of the 
deposition rate differences in the Ar-lg and He-lg (g being the terrestrial 
gravitational acceleration) cases need to be done with respect to the baseline 
(Og) expected difference (^13% for Tsui's experiments). Furthermore, this 
approach will not be able to detect the subtle self-compensating effects of 
buoyancy resulting from complex couplings, i.e. even when the difference 
between the deposition rates for the terrestrial Ar and He experiments is the 
same as the ones expected in the absence of gravity, one can still not 
conclude that the effect of buoyancy on the deposition process is negligable. 
This study, while assessing the type and the strength of convection in both 
terrestrial cases, also compares the lg cases to the respective ideal Og 
cases, and, hence, eliminates the ambigiuties in possible interpretations. 
This, of course, is one of the advantages of a computational approach. 

Because of the steep temperature gradients (V700K/cm) and complex reactor 
geometry, simple dimensionless analysis using quantities such as the Grashof 
and the Rayleigh numbers fail to provide reliable information for the degree 
of significance and the specific effects of buoyancy. However, they can still 
be used as a guide to gain some insight for the expected magnitude of free- 
convection in various cases. Consequently, the Ar-lg case is expected to have 
larger buoyancy effects mostly because it is heavier than He and the buoyancy 
effects for the He-lg case is expected to be significantly suppressed. A 
comparison of the fully calculated flow fields for the Ar-lg and the Ar-Og 
cases, performed by examining Lagrangian particle trajectories in Figure 3, 
shows the dramatic 3-D influence of buoyancy. For this purpose neutral- 
buoyant and non-inertial particles are introduced into the stream at 0.3cm 
above the susceptor level and their trajectories are calculated. Such 
numerical "flow visualization" tools clearly demonstrate the "cork-screw-like" 
behavior of the flow which the 2-D analyses cannot be expected to correctly 
describe. It should be noted that the three dimensionality of the flow inside 
the reactor can be observed even for the Ar-Og case due to the cross-sectional 
flow-area changes, volumetric expansion of the gases due to thermal heating, 
etc.. On the other hand, a close comparison of the He-lg and He-Og cases in 
the present study indicates that the buoyancy effects are, indeed, 
significantly suppressed. It is noted that such a change is accomplished by 
simply replacing 25 mole percent Ar with 25 mole percent He in a 75 mole 
percent H 2 mixture. The subtle differences between the He-lg and He-Og 
cases and how they may be important in some CVD applications are beyond the 
scope of the objectives of this paper and will be elaborated upon elsewhere 
(Gokoglu et al., to be published). Therefore, we will omit the presentation 
of the He-Og cases in this specific study and point out, again, that the 
overall flow features in both cases result in similar transport phenomena, and 
consequently, similar Si deposition rates. 
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Figure 4 shows the velocity vectors at the previously designated three 
X-slices for the three cases: Ar-lg, Ar-Og and He-lg. The complex flow 
pattern for Ar-lg should be contrasted with the other two and, especially, 
with the smooth behavior of Ar-Og where there are no vectors pointing opposite 
to the main flow direction (reverse flow). The comparison of Ar-lg and He-lg 
demonstrates the intensity of buoyancy-driven convection rolls for Ar-lg and 
how it influences the flow on the susceptor. In fact., it has been 
computationally determined, although not readily observable from Figure 4, 
that the number of f ree-convection cells across the width of the reactor 
varies along the susceptor length. The presence of the vertical drops on both 
sides of the susceptor helps to stabilize the flow on the susceptor, as is 
evident from the He-lg case, by localizing the gravitationally-induced rolls 
to the drop region. The He-Og case (not shown) is qualitatively like Ar-Og 
with no reverse flows or f ree-convection cells and, therefore, is 
fundamentally different from the He-lg case. However, as is mentioned above, 
buoyancy effects are so substantially supressed for He-lg in this study that 
the presentation of He-Og results is found to be redundant with .respect to the 
He-lg case. 

Figure 5 shows Z-slices of the modelled reactor volume for velocity, 
temperature and silane concentration fields. Except for cases (b) and (c), 
these slices are taken 0.19cm away from the center symmetry-plane. For cases 
(b) and (c), however, they are taken 2.10cm away from the center symmetry- 
plane and are slightly beyond the susceptor. The cases (a), (b) and (c) 
depict the existance of a convection cell also in the longitudinal direction, 
which is stronger for Ar-lg than for He-lg. Case (a) shows that the reversed 
flow in the lower portion of the reactor exit zone rotates towards the top of 
the reactor and significantly affects the axial main flow at the trailing edge 
of the susceptor. Case (b) shows that the reversed flow extends further 
upstream along the sides of the susceptor, turn up and intensify due to the 
heated susceptor side walls, and, eventually, turn around towards the reactor 
exit. It is interesting to observe that while the flow naturally "falls” 
after every geometric down-step for Og cases, as one would intuitively expect, 
it actually rises for Ig cases. The presence of these longitudinal cells on 
both sides of the susceptor implies, especially for Ar-lg, that the main axial 
flow on the susceptor is likely squeezed from both sides and accelerates 
through the reduced effective flow-area. Case (c) demonstrates the apparent 
supression of this buoyancy related phenomenon by helium. 

Cases (d), (e) and (f) in Figure 5 are isotherm contour plots for Ar-Og, Ar-lg 
and He-lg displaying also the calculated heating of the surface upstream of 
the susceptor and the temperature distribution inside the solid bodies. Note 
that the susceptor is kept at a constant temperature of 1373K. The comparison 
of (d) and (e) shows how the temperature field is modified by the f ree- 
convection effects, whereas the comparison of (e) and (f) shows the relative 
effect of different Prandtl numbers (the ratio of the kinematic viscosity, 

Vi/p, to the thermal dif f usivity ,' A/(pc p ), of the gas mixture) on the 
temperature distribution. 

Cases (g), (h) and (i) display the relative effects of convection and 
diffusion on silane depletion (related to Si deposition) rates. The presence 
of f ree-convection facilitates the convective mass transport to the deposition 
surface resulting in more depletion of the source gas silane. However, 
because the diffusional mass transport is very significant for the delivery of 
Si-carrier species to the surface in this study, using a He host gas mixture 
depletes more silane even with supressed f ree-convection. 

Figures 6, 7 and 8 are 2-D contour plots of temperature and silane, silylene 
and disilane concentrations around the susceptor at the three X-sliced 
locations along the susceptor for the Ar-lg, Ar-Og and He-lg cases, 
respectively. The isotherms in Figure 6 demonstrate clearly how free- 
convection distorts the temperature profiles around the susceptor for Ar-lg. 

As the thermal field develops and the number of f ree-convection cells varies 
over the susceptor along the flow axis, the nature of the distortion also 
varies, both qualitatively and quantitatively. In fact, for Ar-lg the thermal 

ORIGINAL PAGE IS 
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field never becomes fully developed over the susceptor in this particular 
study , although for Ar-Og and He-lg it does. The consequent appearance of the 
species concentration fields is influenced in two major ways due to these 
buoyancy effects. Firstly, the rates of homogeneous reactions that produce 
the reactive species respond exponentially to the local temperature; and, 
secondly., the free-convective transport of species both in the Y-Z planes 
( cross-sectional cells) and in the X-Y planes (longitudinal cells) affects 
both the local homogeneous reaction and transport rates. Note that when the 
buoyancy effects are eliminated or supressed (Figures 7 and 8) the problem can 
be analyzed fairly accurately by (pseudo) 2-D treatments which consider the 
flow and thermal field developments and their subsequent effects on the 
concentration fields using boundary layer approximations. In such cases, 
cross-sectional and longitudinal mixing will depend basically on diffusion. 

A comparison of SiH 4 concentration contours among different X-slices in 
Figures 6(b), 7(b) and 8(b) indicates that the concentration of silane gets, 
smaller closer to the hot substrate surface. This is attributed to the 
formation and faster surface reaction (deposition) of the other reactive 
species, SiH 2 and S^Hg, rather than to the relatively much slower surface 
reaction (deposition) of SiH^ at the susceptor temperature of this study. 
Indeed, this assertion is confirmed by the further analysis of numerical data 
which indicates that, among the three Si-carrier species, the contribution of 
the SiH 4 species to the total Si deposition rate is negligably small as 
compared to the other two. As can be seen from the comparison of Figures 6(b) 
and 7(b), the depletion of silane is more for Ar-lg than for Ar-Og due to the 
enhanced transport of the source to the reaction zones by f ree-convection. In 
fact, the uniform cross-sectional silane concentration observed at the reactor 
exit-plane for Ig cases is the result of such more efficient mixing, although 
a silane concentration gradient still exists across the reactor exit-plane for 
Og cases. 

The dissociation reactions for SiH 4 and Si 2 H^ are favored at higher 
temperatures. Therefore, the formation of S 1 H 2 is promoted closer to the 
hot substrate surface. It is observed, although not readily seen from Figures 
6(c), 7(c) and 8(c) due to the insufficient resolution, that SiH 2 has a 
concentration peak close to the hot surface. It gets less concentrated in the 
cooler regions due to slower SiH 4 dissociation and faster formation, 

and less concentrated in the immediate vicinity of the hot surface due to the 
instantaneous surface reaction, resulting in positive diffusional flux toward 
the surface. 

The formation of S^Hg , Figures 6(d), 7(d) and 8(d), in cooler regions forlg 
cases is promoted by the f ree-convection facilitated transport of S iH 2 , which 
has been formed in hotter regions. The reader is alerted to the fact that 
S^Hg contour values are one to two orders of magnitude smaller than SiH 2 
values.- Concentration contours for S^Hg at the susceptor leading edge show 
that, similar to SiH 2 , a maximum exists also for S^Hg. However, subsequent 
X-slices downstream of the susceptor show that S^Hg concentration uniformly 
decays toward the susceptor top surface, suggesting that the homogeneous 
reactions reach equilibrium and the S^Hg concentration field is governed by 
the transport rates with an infinately fast surface reaction (relative to 
transport rates). It should be noted that although the gas phase chemical 
reaction mechanisms are the same for both Ar and He cases, the significance of 
gas phase reactions are different for each case depending on the relative 
rates of convective and diffusive transport. 

Figure 9 depicts the contour plots of resulting Si deposition rates obtained 
from our computations in comparison to Tsui 1 s/16/ corresponding experiments. 

It should be mentioned that the experimental contours are generated from a 
total of only twelve data points for illustrative comparison purposes. The 
measurements were made at the intersections of four equally-spaced rows and 
three equally-spaced columns outlined in the central region of the susceptor. 
Therefore, the shapes of the experimental contours should be interpreted 
accordingly. All of the three calculated deposition rate contours show the 
same trend of decreasing rates toward the center and trailing edge of the 
susceptor, resembling an amphitheater. Regardless of the level of free- 
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convection, these shapes reflect a) the depletion of Si-carrier source species 
by deposition toward the trailing edge, resulting in lower rates, and b) the 
steeper gradients of solution variables toward the sides due to the geometric 
and thermal setup, resulting in higher rates. The widthwise variation of 
deposition rates is affected by the cross-sectional f ree-convection cells in 
two ways : 1). the rates tend to be further reduced at the center symmetry- 
plane and enhanced at the sides because the cells rotate up toward the top 
wall of the reactor at the center and down along the sides, and 2) the rates 
tend ,to. be augmented, though unevenly, everywhere along the width by the 
facilitated supply of reactant species from higher concentration regions to 
deprived areas. These two competing effects result in the overall enhancement 
of rates. Figure 9(a) and (c), with varying magnitudes along both the width 
and length of the susceptor. 

A comparison of integrated deposition rates over the total surface of the 
susceptor .for Ar-lg, He-lg and Ar-Og, Figure 9(a)-(c), indicate that the 
highest rate is obtained for He-lg and the lowest rate is for. Ar-Og. 
Considering, that f ree-convection is significantly suppressed for He-lg, this 
shows that the faster diffusion rates of Si-carrier species in He' (versus in 
Ar ) has a more pronounced effect in enhancing rates than f ree-convection. 

A comparison of computational and experimental deposition rate results. Figure 
9(a) and .(b), reveals that our preliminary modelling effort, is capable of 
predicting the qualitative and quantitative behavior of Si deposition from 
silane with a remarkable accuracy. Although the axial trends are more 
correctly predicted, the difference in widthwise variation can partially be 
attributed to the scarcity of experimental data. 


4 - CONCLUSIONS AND CLOSING REMARKS 


The computational fluid dynamics code FLUENT-CVD has been adopted to simulate 
a 3-D experimental CVD reactor. The code incorporates the effects of both 
homogeneous (gas phase) and heterogeneous (surface) chemistry. It has been 
applied to silicon deposition from silane at a substrate temperature of 1373K 
using the finite rate homogeneous reactions relevant to silylene and disilane. 
Either a 25% helium or a 25% argon by mole with hydrogen balance is used as 
the host gas mixture for comparison with experiments. These two cases were 
shown . to have appreciable (Ar mixture) or supressed (He mixture) levels of 
gravitationally-induced f ree-convection. The simulations were also carried 
out in the absence of gravity to eliminate f ree-convection. The complex 3-D 
nature of the flow under the terrestrial gravitational field has been 
demonstrated. The reported results are considered to be preliminary due to 
the current limitations of the code and the relatively crude computational 
mesh used. A more thorough analysis will be published later. However, the 
correct physicoand thermo-chemical phenomena have been sufficiently resolved 
in the present study. The numerical results compare favorably with 
experiments. 

For the case presented here, a major enhancement of the code will be the 
inclusion of thermal (Soret) diffusion. A more meaningful code verification 
study with a finer computational mesh can then be carried out using Tsui* s/16/ 
data for various Ar-H 2 anc ^ He-H 2 mixture ratios at different substrate 
temperatures, in the absence of homogeneous nucleation. The capabilities of 
the code can be assessed for the extreme limits where there are no homogeneous 
and/or heterogeneous chemical reaction barriers. The effects of changes in 
the reactor geometry, orientation with respect to the gravity vector, 
uncertainties in chemical kinetic parameters, etc. will be the subject of our 
future publications . 

The improved design and optimization of future CVD reactors requires improved 
and optimized computer codes. It is essential that the data fed to such 
computational codes be of utmost quality. However, the refinement of such 
input information can more efficiently be achieved via more focused 
experimental and theoretical efforts. Therefore, modelling efforts addressing 
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specific aspects of a complex phenomena in a better defined environment/24/ , 
or smaller scale experimental efforts geared towards providing answers to 
individual thermo-chemical questions related to CVD are integral parts of a 
coordinated program. Only when such a loop is closed can computational codes 
go beyond research purposes providing invaluable insights and be confidently 
used' as powerful design tools. 
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assistance of Mr. Ron Gaug (Univ. of Akron, Akron, Ohio). Special thanks are 
also due to Mr. Kurt Fretz, Drs. Tom Jasinski and Zahed Sheikholeslami 
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LIST OF SYMBOLS 

A pre-exponential factor for Arrhenius reaction rate expression. 

Cp specific heat at constant pressure per unit mass of gas mixture. 

D Fick diffusion coefficient. 

E activation energy for Arrhenius reaction rate expression, 

g earth gravitational acceleration, 

h enthalpy of gas mixture. 

J total number of gas phase species, 

k reaction rate constant. 

K total number of reactions, 

n unit vector normal to surface, 

p pressure. 

R universal gas constant. 

Rjk mass rate of creation of species j per unit volume (gas phase) 
s by reaction k. 

Rjk mass rate of creation of species j per unit area (surface) by 
reaction k. 

Sc Schmidt number, y/(pD). 

T temperature in Kelvin degrees, 

u velocity. 

W molecular weight, 

x dimension coordinate. 

Y mass fraction. 

3 exponent of T in Arrhenius reaction rate expession. 

y surface reaction probability. 

A thermal conductivity of gas mixture. 

M dynamic viscosity of gas mixture, 

p density of gas mixture. 

3 partial differential operator, 

v vector differential operator. 

Subscripts 

i coordinate direction, 

j pertaining to species j . 

k pertaining to reaction k. 

l dummy variable. 

Superscripts 

s pertaining to surface. 
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TABLE I : Contour values defined in Figs. 5-8 


MASS FRACTION 


Contour 

Temp . , K 

SiH 4 

SiH 2 

Si 2 H 6 

A 

1 . 33 ( +3 ) 

1 - 2 ( — 2 ) 

3 . 4 ( -4 ) 

2 . 0 ( -5 ) 

B 

1 . 06 ( + 3 ) 

9 . 0 ( -3 ) 

2 . 8 ( -4 ) 

l.M-5) 

c 

7 . 92 ( +2 ) 

5 . 8 ( -3 ) 

2 . 2 ( -4 ) 

4 . 0 ( -6 ) 

D 

5.24(4-2) 

2 . 7 ( -3 ) 

1 . 6 ( -4 ) 

3.1(-6) 

E 

— 

1 . 8 ( -3 K 

1 . 4 ( -4 ) 

1 . 0 ( -6 ) 

F 

— 

l.M-3) 

1 • 0 ( -4 ) 

2.0(-7) 

G 

— 

3 - 5 ( -4 ) 

2 . 0 ( -5 ) 

— 

H 

— 

— 

1 . 5 ( -5 ) 

-- 
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Figure 1. 

The schematic of the modelled portion 
of the reactor. 



Figure 2. 

Magnified view of the susceptor section, 
showing three x-slices at 0.0cm, 3.58cm 
and 7.0cm from the susceptor leading edge. 
Data in Figures 4 and 6-9 are presented 
according to this format. 


ml 



Figure 3. 

Trajectories of particles introduced 0.3 cm above the surface at 
the entrance of the modelled reactor volume for (a) Ar-Og, and 
(b) Ar-lg. 
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Position from 
center symmetry 




















Figure 7. 

The x-slices for Ar-Og: (a) isotherms, and concentration contours 
for (b) SiH^ (c) SiH 2 and (d) S^Hg. Refer to Table I for 
contour values . 
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